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We formulate a smoothed-particle hydrodynamics numerical method, traditionally used for the 
Euler equations for fluid dynamics in the context of astrophysical simulations, to solve the non-linear 
Schrodinger equation in the Madelung formulation. The probability density of the wavefunction is 
discretized into moving particles, whose properties are smoothed by a kernel function. The tradi¬ 
tional fluid pressure is replaced by a quantum pressure tensor, for which a novel, robust discretization 
is found. We demonstrate our numerical method on a variety of numerical test problems involving 
the simple harmonic oscillator, soliton-soliton collision, Bose-Einstein condensates, collapsing sin¬ 
gularities, and dark matter halos governed by the Gross-Pitaevskii-Poisson equation. Our method 
is conservative, applicable to unbounded domains, and is automatically adaptive in its resolution, 
making it well suited to study problems with collapsing solutions. 
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I. INTRODUCTION 

Quantum mechanics is one of the basic pillars of mod¬ 
ern physics. The Schrodinger equation describes the 
quantum mechanical evolution of the wavefunction of a 
particle over time. The non-linear Schrodinger equation 
(NLSE), also called the Gross-Pitaevskii equation, is a 
non-linear extension of the Schrodinger equation, which 
describes the ground state of a quantum system of iden¬ 
tical bosons using a single-particle wavefunction approx¬ 
imation and a pseudopotential model for interaction. It 
is ideal for describing a Bose-Einstein condensate (BEG): 
dilute gas of bosons in a low-temperature state very close 
to absolute zero. BECs were first predicted in the early 
days of quantum theory by Bose and Einstein in 1924- 
1925. The first realization in the laboratory was achieved 
in 1995 [1, 2], which marked a new era in atomic, molec¬ 
ular and optical (AMO) physics and quantum optics [3]. 
The NLSE has applications and extensions to entirely 
different physical systems as well, including the propaga¬ 
tion of light in non-linear fiber optics [4] , Langmuir waves 
in plasmas [5] , and self-gravitating BEG models for dark 
matter, governed by the Gross-Pitaevskii-Poisson equa¬ 
tions [6]. 

The NLSE is challenging to solve and almost always 
requires numerical solutions. Ongoing research has led 
to the development of a variety of methods to solve these 
systems in various contexts, such as those for solving 
time-evolution of BEG systems [3, 7-10] and obtaining 
their ground states [11-14]. These methods solve for the 
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solution to the NSLE in the standard form, and typi¬ 
cally employ finite-difference, finite-element, or spectral 
methods. Other non-standard methods for solving quan¬ 
tum systems have been proposed as well, such as lattice 
Boltzmann [15, 16] and unitary qubit lattice algorithms 
[17, 18]. Each method has different strengths and limi¬ 
tations when applied to different systems [19]. 

We propose a novel, conservative numerical approach 
for solving the NLSE that is quite different from the stan¬ 
dard approaches. We solve the NLSE in Madelung hydro- 
dynamic form, using a smoothed particle hydrodynam¬ 
ics (SPH) algorithm. The Schrodinger equation, as well 
as the NLSE, can be reformulated under the Madelung 
transformation to take a different form that resembles the 
fluid equations [20]. The equation in Madelung form de¬ 
scribes the evolution of the quantum probability density 
of the wavefunction under a quantum “pressure” tensor, 
and is equivalent to the standard form. 

SPH is a particle-based method for computational fluid 
dynamics. It was originally invented to simulate poly¬ 
tropic stellar models under non-axisymmetric conditions 
[21, 22]. It has since been extended and coupled with 
additional physical processes and plays a central role 
in astrophysical and cosmological simulations [23-25]. 
SPH operates independently of any grid, unlike finite- 
difference, finite-volume, or hnite-element methods, and 
interactions between volume elements, such as the pres¬ 
sure gradient, are represented as a force between par¬ 
ticles. The method is purely Lagrangian, meaning that 
interactions and derivatives are evaluated in a coordinate 
system attached to a moving fluid element. The two fun¬ 
damental ideas of SPH are (1) to evolve the positions 
and velocities of particles according to the calculation of 
the forces on each particle at each time step, and (2) to 
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use an interpolating/smoothing kernel to calculate forces 
and spatial derivatives. 

SPH has some desirable inherent features for quantum 
systems. The method is conservative, so the normal¬ 
ization condition on the wave function is preserved to 
machine precision. Also, the SPH method also has no 
domain restrictions; the wave function is free to travel 
anywhere in physical space, which is an advantage grid- 
based methods do not possess. The Lagrangian nature 
of SPH also makes the method useful for the study of 
highly dynamic solutions, such as rapidly rotating wave- 
functions. Furthermore, the SPH method is automati¬ 
cally adaptive in its resolution, making it possible to eas¬ 
ily resolve collapsing features in the solution, and this is 
one of the primary reasons it often is used in cosmological 
simulations of structure formation. 

Our method may be extended in a relatively straight¬ 
forward manner to various other quantum systems, such 
as multi-component, rotational, dipolar, or spin-orbit 
coupled BECs; to Euler-Korteweg systems for capillary 
fluids; or to the study of the semi-classical limit of the 
Schrbdinger equation [26, 27]. Such applications are left 
to future study. 

The remainder of this paper is organized as follows. 
In Section H we discuss the theory of the NLSE and the 
Madelung formulation. In Section HI we describe our 
numerical SPH method for the NLSE. In Section IV we 
demonstrate the accuracy and results of our method on 
some numerical test problems. In Section VI we offer 
concluding remarks and list physics areas of applications 
for the numerical method. 


II. THEORY 


where p = is the quantum probability density of the 
wavefunction, and u = V0, where ijj = |'i/'|exp (id(x, t)). 
The variable 

P = —ipV®Vlnp (5) 

is the quantum pressure tensor. Often in the literate, the 
quantum pressure term is instead written in terms of a 
quantum potential Q = but it turns out that 

for the purposes of discretizing the equations to obtain 
an SPH method, the pressure tensor formulation is more 
useful. 

Optionally, one may add artificial damping to the 
equations, with damping parameter 7 , as 

9tu-1-u • Vu = —-VP —-W — yu. (6) 

P P 2 

This can be useful to bring solutions to a steady state. 
Such a term describes dissipative quantum systems [20], 
where the system loses energy with time. For our pur¬ 
poses, the damping term is useful in relaxing an arbitrary 
wave function to it’s ground state. 

The NLSE has been coupled with self-gravity to form 
the Gross-Pitaevskii-Poisson equation. These equations 
describe models for BEG dark matter halos [ 6 ]. In this 
case, the potential is computed from the wavefunction 
according to Poisson’s equation; 

VV = M^inGp. (7) 

(Note, there is a factor of M^, where M is the total mass 
of the system, because we are using units where p = 
is dimensionless.) 

Another variant of the NLSE is the focusing NLSE 


The Schrodinger equation for quantum mechanics may 
be written in dimensionless form {h = 1 ) as 


idttp{x, t) 


2 


V(x) 


ijj, X e 


t>0. ( 1 ) 


The dynamics of a BEG is well described by the NLSE. 
The NLSE in dimensionless form [ 8 ] may be written as 


idti> = 


2 


V + gW 


i’- 


( 2 ) 


t) = [-V" - X G (8) 

where solutions exist that self focus and become singular 
in finite time for the case ad > 2. Numerically solving 
such a blowup solution is challenging because the spatial 
and temporal gradients grow arbitrarily large while small 
perturbations may arrest the critical collapse. Standard 
grid methods break down and more sophisticated meth¬ 
ods, which involve dynamical rescaling, have been re¬ 
sorted to in order to handle blowup solutions [28]. 


In this form, the normalization is /j'0j^d^x = 1. g G K 
is treated as an arbitrary dimensionless parameter which 
measures the strength of nonlinear interactions. 

Under the Madelung transformation [20], the NLSE 
resembles the fluid equations: 

dtp + y- (pu) = 0, (3) 

atU-Hu-Vu=--V • P--V—- VV, (4) 

p p 2 


III. NUMERICAL METHOD 

We discretize the quantum probability density p = 
of the wavefunction as a collection of N particles. Each 
particle is assigned a “mass” ruj = 1/N so that we satisfy 
the normalization condition ~ = 1 to 

machine precision. 

We wish to solve for the dynamics of p using the 
Madelung formulation. Equation 3 is just a statement 
about the conservation of the normalization condition, 
which is automatically satisfied by our discretization. 
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Equation 4 describes the equation of motion of the par¬ 
ticles. The left-hand side is a convective (Lagrangian) 
derivative of the velocity: ^, namely, the acceleration. 

The main goal of the SPH method is to evaluate the 
acceleration of each particle and update the velocities 
and positions of the particles with each time step using 
an integrator scheme (Section HID). 

With SPH, the value of a field at any point in the do¬ 
main is obtained by smoothing out the values associated 
with the particles. Consider the (trivial) identity: 


A. Calculating density 

Given a distribution of N particles in physical space, 
the density pi at each particle i is required to esti¬ 
mate any field quantity, because it appears in Equa¬ 
tion 14. This is calculated straightforwardly by substi¬ 
tuting H(x) = p(x) into Equation 14: 

p^ = Y^mjW^j (16) 

3 


A(x) = J A(x')5(x — x') dx' (9) 

where A(x) : i—>• K. is any arbitrary function and 

(5(x) is the Dirac delta function. In the SPH scheme, 
S is replaced with an approximation: a smoothing kernel 
IE(x; h), where h is the smoothing length scale. The 
smoothing kernel must have the properties 


J IE(x; h) d^x = 1, 

(10) 

lim IT(x; h) —(5(x). 

h—iO 

(11) 


and must be non-negative and invariant under parity. 
For our purposes, we choose the Gaussian kernel: 

W(x; h)=(^-j^ye^p{-\\^r/h^) (12) 

where /i is a smoothing-length parameter. Alternate 
choices include cubic-splines, which have compact sup¬ 
port and can make pairwise-interaction calculations more 
efficient. For our purposes, we will most often use a fixed 
value for h, but more sophisticated formulations of SPH 
exist which use adaptive values based on particle number 
density [24, 29] (see Section HIG). 

Hence, an approximation to the field A(x) in Equa¬ 
tion 9 is 

A(x) ~ y A(x')IE(x - x') dx' (13) 

To this equation, we apply a second approximation, 
namely discretization: we sum over the N particles. The 
SPH approximation to A(x) is thus 

where Xj is the location of particle j in physical space 
and pj = p{x-j) (calculated via Equation 16). 

Similarly, gradients of fields can be approximated as 
follows 

VA(x) ~ ^A{xj)VW{x - Xj- h). (15) 

, Pj 


where we have defined 

Wi, = W(x,-x,; h). (17) 


B. Calculating pressure tensor 

In classical fluid dynamics applications, the pressure 
Pi at a particle location is calculated using an equation 
of state, which depends on quantities such as the density 
and/or the internal energy of the fluid. One example is 
the polytropic equation of state P = where k is 

a constant and n is the polytropic index. 

In the case of quantum mechanics, the pressure is 
replaced by a symmetric pressure tensor (Equation 5), 
which is a non-local quantity because it depends on gra¬ 
dients of the density. 

First, the derivatives of the density field are calculated 
at the location of each particle, along the lines of Equa¬ 
tion 15: 


dxPi = Y . (18) 

3 


Similarly for dz- 

The second derivatives can be calculated in similar 
fashion as 


P^xyPi — ^ ^ '^33jjdxy^^ij , 

3 

and similarly for dxx, dxz, dyy, dyz, dzz- But this is not 
the only discretization. They may also be calculated us¬ 
ing the (more accurate) second-order discretization [30]: 

^xyPi — ^ ^ " (Pj Pi) dxyAPij ^ (^0) 

, P3 


Other alternate discretizations exist in the literature as 
well, such as the difference scheme, and are widely used 
for applications such as heat conduction [31, 32]. 

Finally, the components of the quantum pressure ten¬ 
sor of Equation 5 are computed as: 


p — 

' i,xy — 


E 


^1 r {dxPj){dyPj) 
Pj 4 [ Pj 


dxyPj 




( 21 ) 


Note that the gradient operator shifts to the kernel, 
whose derivative is analytically known. 


and similarly for the rest of the components of the tensor. 
This is on of the main equations of our paper, which 
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provides a robust discretization for the quantum pressure 
tensor that yielded well-behaved solutions in all our test 
problems. 

We note here an analogy between the quantum pres¬ 
sure tensor P = —\pS/ 0 Vlnp and the equations for 
inviscid capillary fluids. Such classical, non-ideal fluids 
can be described by the Euler-Korteweg equations, which 
have the form: 


Optionally, one may artificially damp the solutions by 
adding the damping term to the right hand side of 
Equation 24. Damping is useful for obtaining steady- 
state solutions: eigenstates and/or ground-states of the 
system. Is is also useful for generating initial conditions. 

D. Leap frog time integration 


^ = V [n{p)^^p + \k'{p)\S/p\^'^ (22) 

Taking k{p) = l/(4p) leads back to the NLSE. In the clas¬ 
sical case, the pressure tensor would act as a surface ten¬ 
sion term, and its sign would act to add decoherence into 
the solution (which is a typical quantum phenomenon). 
An opposite sign would lead to cohesion. 


C. Calculating acceleration 


A key guiding principal in formulating an SPH method 
for obtaining robust results is to choose discretizations for 
the forces the particles experience such that the forces 
between pairwise particles obey Newton’s third law, i.e., 
are equal and opposite. This allows for the particles to 
quasi-regularize as they sample the true solution of the 
field [25]. This leads to better-than-random-Monte-Carlo 
sampling in the reconstruction of field quantities. 

In the standard SPH formulation with scalar fluid pres¬ 
sure, the acceleration of a particle due to the pressure 
gradient — is calculated in symmetric fashion as: 


du., 

dt 


j 




(23) 


The particles are initialized with positions and veloci¬ 
ties dictated by the initial conditions of the problem (in 
problems that use damping to reach a steady state, the 
initial conditions can be chosen to be a random or uni¬ 
form distribution of particles with 0 velocity; the particle 
configurations of such steady state solutions can then be 
used for problems with dynamics). Then, equation (4) 
may be solved with a time integration method, such as 
Runge-Kutta or leap frog. The leap frog method is often 
preferred because it is explicit and symplectic. In our im¬ 
plementation, we use the second-order leap frog scheme 
as follows: 


u{t -I- At/2) = u(t — At/2) -I- a{t)At (25) 

x(t + At) = x(t) -f- u(t + At/2)At (26) 

that is, we calculate positions and velocities at inter¬ 
leaved time points. We note that at the start of the 
simulation we only know the initial conditions x(0) and 
u(0), and must use first-order Euler to step back half a 
time step and find u(—At/2). 

To find the velocities at the same time step intervals 
as the positions, we use the approximation: 

u(t) = i (u(t — At/2) -I- u(t -I- At/2)). (27) 


The force due to the pressure tensor can also be calcu¬ 
lated in a similar fashion: 


dui 

dt 


= -E’ 


[^i,xx! Pi,3:2] ^j,xyi PjjXz] 


^2 

■ + 

^2 

Pi 


Pj 

[^i,yx, ^i,yyi ^i,yz] 

■ -1- 

[Pj,i/X) Pj,yy; ^j,yz\ 

Pi 


[Pi,22:; Pi,2y; Pi,22] 

+ 

[Pj,zx, Pj,zy, Pj.zz] 


Pt 


(24) 


Obeying Newton’s third law between pairwise particles 
is one of the key reasons we chose to discretize a version 
of the Madelung equations that uses the pressure tensor. 
Note the inner square brackets here denote a vector, and 
we have written out the entries explicitly. 

We can also calculate the additional acceleration due 
to the non-linear term — in Equation 4 by setting 


P, = 


gp 


2 in right-hand side expression of Equation 23. 


E. Pseudocode 


A pseudocode of the main loop of the SPH algorithm 



is shown below. 

■yWij, 

Main Loop 

S7Wij, 

for t = I : N_time_step 


% leap frog 

yWij, 

v_phalf = v_mhalf -t- 
x -|-= v_phalf * dt; 


V = 0.5 * (v_mhalf 4- v_phalf); 
v_mhalf = v_phalf; 

% update densities, pressures, accelerations 
rho = GetDensity(x, m, h); 

P = GetPressure(x, rho, m, h); 

a = GetAcceleration(x, v, m, rho, P, b, beta, lambda, h); 

end 

Our implementation is 0{N‘^) because we calculate all 
pairwise interactions between particles. Alternative im¬ 
plementations have been developed in literature to make 
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the computations more efficient {0{N log N).), such as 
computing interactions only between the k nearest neigh¬ 
bors with using a hierarchical tree algorithm [33]. 


F. Gross-Pitaevskii-Poisson 

The resulting equations of motion from self-gravity 
(Equation 7) can be calculated using and 7V-body tech¬ 
nique (here M = 1): 


^ ^ ^ rrij jrjr,) 

dt ^ (|r,r,|2 +£2)3/2 1 ) 

where e is a smoothing length, used to avoid numerical 
problems of close particle encounters (where the acceler¬ 
ation blows up) in collision-less dynamics. In SPH, e is 
typically equated with the kernel’s smoothing length h. 
SPH couples naturally with the Poisson equation, mak¬ 
ing our method well-suited for finding solutions of the 
Gross-Pitaevskii-Poisson equations. 


G. Adaptive Smoothing Lengths 


In some applications it is advantageous to use an adap¬ 
tive smoothing length hi for each particle i for improved 
numerical accuracy. The variable smoothing length de¬ 
pends on the density at the fluid and allows the algorithm 
to handle regions with high and low densities more pre¬ 
cisely [34]. 

The adaptive smoothing length and the density can be 
calculated self-consistently using an iterative scheme as 
follows. The density estimator becomes: 

Pi = ^ mjW{yi.i - xj; hi) (29) 

j 

There is a correction factor to the momentum equation 
to allow for the spatial variation in smoothing lengths in 
the equation of motion (Equation 23) 


dt 



Pj 


VIE(xi - Xj] hi). 


(30) 


where (for 3D case) 



dW (x,; — Xj ; 

dhi 


hi) 


(31) 


An analogous correction correction applies to Equa¬ 
tion 24. 

The value for hi is determined by solving for the root 
of 


C{hi) = ruj 



hi) 


(32) 


using a Newton-Raphson iterator (or alternate tech¬ 
nique). p is an order unity constant; for our applica¬ 
tions, we use rj = 1.4. This constraint ensures that 
Pih^ = const, i.e., it maintains a constant mass within 
the smoothing kernel. Given an initial guess for /i^, the 
Newton-Raphson iterator gives an updated value hi^new 
according to: 



until a tolerance threshold is reached: J/ii^new — hi\/hi < 
etoi- We choose etoi = 10“^. 


IV. RESULTS 


Here we present a number of simple numerical tests 
to demonstrate our SPH method for the Schrodinger 
equation and the NLSE. The aim of these tests is to 
demonstrate the accuracy of our method (in capturing 
steady-state solutions and dynamics) and highlight differ¬ 
ent physical regimes well-suited for our numerical method 
to handle (such as systems with self-gravity or collapsing 
solutions). 



N 


FIG. 1. Convergence properties of the SPH code for obtaining 
the ground state of the Schrodinger equation with a ID sim¬ 
ple harmonic oscillator potential. Second order convergence 
is achieved by increasing particle number N (and decreasing 
smoothing length as l/N). The obtained ground state solu¬ 
tion for p at the various resolutions is shown in the inset (the 
shades of teal correspond to the shades of the circles in the 
convergence plot: higher resolution approaches the analytic 
answer, which is shown with the thick gray line). 
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A. ID simple harmonic oscillator ground state 

First, we demonstrate the algorithm’s ability to recover 
the ground state of the Schrddinger equation with a ID 
simple harmonic oscillator potential. This is a very sim¬ 
ple test, but can be used to demonstrate our code’s con¬ 
vergence properties. The particles are initially drawn 
from a uniform distribution in the range [—4,4], and we 
add a damping term with 7 = 4 to relax the system to 
the ground state. We evolve the system in the potential 
V{x) = until a steady-state configuration is found. 
We use a smoothing length oih = 200/A and a time step 
of 0.004. The exact solution is given by 

Po = TT-^/^exp (-x^) (34) 

Figure 1 demonstrates the convergence rate and solu¬ 
tion. This method of obtaining the ground state is also 
useful for generating initial conditions for future tests. 



FIG. 2. (Color online) Time evolution of oscillating wave- 
function in simple harmonic oscillator potential. The SPH 
approach captures the dynamics well. The analytic solution 
is shown in thin gray lines. 


B. ID simple harmonic oscillator dynamics 

We consider the time evolution of a wavefunction in 
the potential V{x) = ^x^ with initial conditions: 

Po = 7 r“^/^exp (-a;^) , (35) 

and 

(36) 



X 


FIG. 3. (Color online) Evolution of two solitons in simple har¬ 
monic oscillator potential. The SPH code captures the profile 
well, which consists of two solitons oscillating in opposite di¬ 
rections. 


The evolution has analytic solution 

p(t) = 7 r“^^^exp (—(a; — sin(t))^) . (37) 


Our simulation used N = 300 particles, smoothing 
length h = 0.2667, and time step dt = 0.01. The ini¬ 
tial conditions are obtained by damping a random con¬ 
figuration of particles to the ground state, as in Sec¬ 
tion IV A. Figure 2 shows the evolution of the wavefunc¬ 
tion with time, which matches the analytic solution very 
well, showing exact agreement with the periodicity of the 
solution. We note the tiny offset in the peak of the wave 
function compared to the analytic result is due to the 
finite number of particles we use, and this discrepancy is 
reduced with increasing particle number and decreasing 
smoothing length. 

We also demonstrate the dynamics of two solitons in a 
harmonic potential evolved under the linear Schrddinger 
equation. Two copies of the initial condition in the pre¬ 
vious setup are superimposed on top of each other. One 
soliton is given initial velocity u = - 1-1 and the other 
V = —1. In this simple example, the solitons do not in¬ 
teract and just pass through each other. This example 
demonstrates the algorithm’s ability to capture multiple 
phase dynamics. Our simulation used N = 800 particles, 
smoothing length h = 1, and time step dt = 0.001. Fig¬ 
ure 3 shows the evolution of the wavefunction with time, 
which matches the analytic solution very well. 


uo = 1 . 

















lytic profile 



FIG. 4. (Color online) Collision of two bright solitons evolved 
under the NLSE. Collision occurs at t = 5, during which 
fringes are formed. The solitons return to their original profile 
after interaction. 



FIG. 5. (Color online) Collision of two dark solitons evolved 
under the NLSE. Collision occurs at t = 2.5. 


C. ID NLSE soliton-soliton collision 


We simulate the collision of two solitons that are so¬ 
lution to the NLSE to show the stability and accuracy 
of the pre- and post- collision nonlinear states with our 
method. We initialize two bright solitons [35] with ana- 


'i/'o(a;) = ^sech(a; ± xo)e^™^ (38) 

(such an initial condition is initialized by damping an ini¬ 
tially uniform distribution of particles under the appro¬ 
priate potential that gives rise to the profile as the steady- 
state solution.) A single bright soliton moves at constant 
velocity w, with peak initially determined hy x = xq. We 
simulate the interaction of two solitons initially located 
at a; = ±a;o = ±5 and each traveling with speed r; = 1 to¬ 
wards each other, evolved under the NLSE with g = —1. 
During interaction, oscillations are created, and at col¬ 
lision the peak value of the amplitude doubles from the 
initial peak values of the single soliton, as expected in 
soliton-soliton collisions [35]. After interaction, the soli¬ 
tons return to their original profile shape and continue 
traveling at velocity v. The results of our simulation, 
demonstrating the mentioned behavior, is shown in Eig- 
ure 4. The simulations use N = 100 particles, smoothing 
length = 1, and time step dt = 0.001. The pre- and 
post collision profile is preserved very well under evo¬ 
lution: our initial condition has profile peak 0.2529 and 
velocity 1, and the post-collision profile at f = 8 has peak 
0.2526 and velocity 1.0043 (note that the analytic solu¬ 
tion has initial peak 0.2500, so there is a small numerical 
offset in the initial conditions due to the truncation errors 
associated with relaxing the SPH particles in a potential 
that yields the initial conditions). 

We point out that, as is inherent in SPH, in regions 
where density is low, the profile may be dominated by 
the shape of the kernel (e.g. the slight bump in density 
at a; = 0, f = 8). Additionally, there are truncation 
errors in the inital conditions generated by relaxing the 
SPH particles in a constructed potential to generate the 
inital conditions. This can potentially lead to additional 
small amplitude waves propagating in the time evolving 
solution. 

We also demonstrate that the SPH method can cap¬ 
ture collision of dark solitons. In the case, the profile we 
consider is 


ipo{x) = iv ± tanh(x ± xq) (39) 

with initial positions determined hy x = ixo = ±5, and 
velocity v = 2. Like the bright soliton, the profile for a 
single dark soliton remains unchanged under evolution of 
the NLSE. The velocity of the SPH particles making up 
the soliton is 


u{x) = T 


xsech^(x ± xo) 

+ tanh^(x ± Xq) 


(40) 


Thus, unlike the bright soliton, the initial Lagrangian 
velocities of the SPH are not constant (note, the soliton 
profile still moves at constant velocity). Thus this test 
demonstrates dynamics and ability to maintain a soliton 
profile under a continuum of phase velocities. The system 
is evolved by the NLSE with g = 1. The results of our 


















simulation are shown in Figure 5. The simulation uses 
N = 100 particles, smoothing length h = 0.4, and time 
step dt = 0.001, and periodic boundary conditions on the 
domain [—10,10]. The simulation preserved the profile 
well after interaction. 

We point out that a general issue with the fluid-like ap¬ 
proach of the Madelung equations is that the quantum 
pressure term can be a problem from the singularity point 
of view. This is certainly an issue for grid-based meth¬ 
ods that attempt to solve the Madelung equations. The 
density can approach zero and the corresponding velocity 
can be infinite (so that the momentum is finite). How¬ 
ever, with an SPH approach the situation is improved. 
The velocities are always calculated at the locations of 
the particles, at which there is always a minimum den¬ 
sity determined by the smoothing kernel. Regions with 
tiny or 0 valued wavefunction have no particles at the 
location representing the solution. 

D. 2D BEC ground states 

We calculate the ground state of a 2D BEC by relax¬ 
ing a random initial condition evolved under the NLSE, 
with damping damping 7 = 4. We calculate the states 
for BECs with g = 0, 10, 50,100, 250, 500 in a potential 
B(x) = i +y^)- The simulations use N = 100 par¬ 
ticles, smoothing length h = 1, and time step dt = 0.1. 

The results are shown in Eigure 6, and are compared to 
the high resolution numerical simulations in [8] , showing 
good agreement. 



FIG. 6. (Color online) Ground state of 2D BEC with various 
values of the parameter g, as obtained by the SPH approach 
and compared to results from high resolution numerical sim¬ 
ulations in [8]. 



t 

FIG. 7. (Color online) Comparison of the blowup solution of 
the focusing NLSE with numerical solutions obtained via our 
SPH method and a second-order finite difference scheme for 
comparison. The SPH solution produces a collapse with the 
right scaling while the finite difference scheme shows focusing- 
defocusing oscillations early on due to discretization effects. 
In the insets, we show the location of the SPH particles at 
three different times, and draw circles around each that have 
radii proportional to the adaptive smoothing length of the 
particle. The inset at the third time frame has an aggregation 
of 12 particles near the singularity. 


E. 2D self-focusing NLSE 

We simulate a blowup solution of the 2D focusing 
NLSE. The initial condition is given by a Gaussian 
'i/i = 2.99e“^^^“''^^^. We use iV = 36 particles and adap¬ 
tive smoothing lengths. Figure 7 shows the evolution 
of the peak of the wavefunction as a function of time, 
which agrees well with the blowup solution scaling of 
the problem |'i/'|max oc (t = 0 corresponds to 

blowup) [28]. In contrast, a finite-difference approach 
shows focusing-defocusing oscillations once the focusing 
becomes too strong [28] (i.e., the oscillatory behavior of 
the peak of the wavefunction shown in Figure 7). Addi¬ 
tionally, such methods require a large number of resolu¬ 
tion elements. Figure 7 also shows the configuration of 
the SPH particles at different times in the simulation, as 
well as their adaptive smoothing lengths. Many of the 
particles cluster at the center because of the collapse, 
which some of the particles are maintained farther out to 
represent the extent of the wave function at these loca¬ 
tions. The inset at showing the particle conhguration at 
the third time frame has an aggregation of 12 particles 
near the singularity. 
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FIG. 8. (Color online) 3D BEC dark matter halo profiles 
calculated with our SPH method for various values of the 
parameter g. The profiles approach the Thomas-Fermi ap¬ 
proximation for increasing values of g. 


F. 3D BEC dark matter halo 

We compute solutions to the Gross-Pitaevskii-Poisson 
equation that describes self-gravitating BEC dark matter 
halos by relaxing an initial condition using a damping 
parameter 7 = 4 . Our simulations use N = 300 particles, 
adaptive smoothing lengths, and units with G = 1. The 
solutions, for various values of g, are show n in Fi gure 8. 
The length scale is normalized by tq = ircjAG in the 
figure. We simulate the cases oi g = 4,10,40,100. In 
the limit of large g, the numerical solution approach the 
Thomas-Fermi approximation: p oc sinc(7rr/ro). 

The Gross-Pitaevskii-Poisson equation describes one 
possible physical model for the non-baryonic dark matter 
that forms a large fraction of the content in our Universe. 
In this model a fundamental scalar field plays the role of 
dark matter, and the model is a competitor to the stan¬ 
dard A cold dark matter model [36]. Large cosmological 
simulations with the BEC model for dark matter have 
been performed on an adaptively-refined mesh to study 
nonlinear cosmic structure formation of gravitationally 
collapsed objects [37, 38]. 

The Gross-Pitaevskii-Poisson equation describes other 
physical systems as well, such as dipolar BECs. 


V. COMPUTATIONAL EFFICIENCY 

Our SPH approach to find solutions to the NLSE main¬ 
tains the simplicity and computational efficiency of the 


original hydrodynamic SPH method. The method only 
requires that the particle positions, velocities, masses, 
and smoothing lengths be stored in memory. More ad¬ 
vanced techniques, standard in the field of SPH, can be 
used to make the method 0{N \ogN), whereas our sim¬ 
ple implementation to calculate pairwise interactions is 
0{N‘^). The SPH technique has successfully been im¬ 
plemented with over 10^*^ particles on standard central 
processing unit (CPU) clusters with the use of Message 
Passing Interface (MPI) routines [39]. In addition, the 
algorithm is well-suited to the modern graphics process¬ 
ing unit (GPU) and GPU cluster architectures [40, 41], 
which have shown an order of magnitude increase in effi¬ 
ciency compared to CPU approaches. These methods can 
simulate a time step of over 10® particles per second. In 
solving the NLSE, the computations per communication, 
are increased due to the computation of a pressure ten¬ 
sor rather than a simple pressure, which boosts the com¬ 
putational efficiency of the original hydrodynamic SPH 
method. A number of numerical methods exist to calcu¬ 
late the self-gravity term, which shows up in the Gross- 
Pitaevskii-Poisson equation, that are 0{N\ogN). These 
include tree and multipole based methods and will make 
the subject of future investigations. 


VI. CONCLUDING REMARKS 

We have demonstrated a new, simple numerical 
method to solve the NSLE using SPH. The method con¬ 
serves the normalization condition on the wavefunction 
to machine-precision. Additionally, the computational 
domain is unlimited, which is very natural for a wave- 
function. The SPH particles that represent the proba¬ 
bility density of the wave function automatically adapt 
to regions where the density is the largest. This makes 
our method ideal for solving collapsing and singular so¬ 
lutions, an area where standard grid methods face diffi¬ 
culties. One limitation of our method is that the the hy¬ 
drodynamic equations and the Gaussian kernels are not 
well-suited for handling systems with singularities, such 
as scalar quantum vortices (at the vortex core the density 
is zero), leading to singular hydrodynamic equations [42]. 
Investigation of such systems is beyond the scope of the 
present work, and may require alternate kernel functions 
to prevent smoothing out the singularities. 

The implementation is relatively simple and easily ex¬ 
tendable to modifications of the NSLE. The numerical 
method can be applied to a variety of physical systems, 
including BECs, nonlinear optics, capillary fluids, dark 
matter that obeys the Gross-Pitaevskii-Poisson equa¬ 
tions, and collapsing singularities. 
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